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Abstract. Interaction between the stellar wind, including 
the solar wind, and the interstellar medium has long been 
the subject of investigation by both astrophysicists and fluid 
dynamicists. This is, first, due to the possibility of comparison 
of physical models for such interaction with the measurements 
performed by Voyager, Pioneer, and Ulysses spacecrafts. On 
the other hand, a complicated structure of the flow containing 
several discontinuities makes it a challenging problem for the 
application of modern numerical methods both in gasdynamic 
and magnetogasdynamic (MHD) cases. 

In the solar wind case, the problem becomes even more 
complicated, since the charge-exchange processes between 
ions and neutral particles must be taken into account. The con- 
tinuum equations are not applicable to the description of the 
neutral particle motion, for their mean free path is much larger 
than the characteristic length scale of the problem. In this case, 
either approximate coupling models or direct Monte-Carlo 
simulation are required. The spatial nonuniformity of the so- 
lar wind and its perturbations and periodicity make the prob- 
lem three-dimensional and nonstationary. From a mechanical 
viewpoint the problem represents the interaction of the uni- 
form interstellar medium and the spherically-symmetric (or 
asymmetric) solar wind flow. We consider various approaches 
used by different authors to solve this problem numerically. 

The presence of the contact surface dividing the two flows 
rises the question of its stability. We discuss the reasons of 
such instabilities and parameters which influence it. 

The presence of the interstellar magnetic field necessitates 
solution of the MHD equations for proper analysis of the 
obtained data. Although the system of governing equations 
remains hyperbolic in this case, the multivariance of the exact 
solution to the MHD Riemann problem makes inefficient its 
application for regular calculations. On the other hand, the 
solution to the linearized Riemann problem is nonunique. We 
discuss the possible ways of applying the Roe-type methods 
and some simplified approaches for numerical solution of the 
ideal MHD equations. One of the difficulties in the solution 
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of the MHD system is the satisfaction of the magnetic field 
divergence-free condition. Different ways to solve this task 
are discussed. If the magnetic field vector in the uniform in- 
terstellar medium flow is not parallel to the velocity vector, the 
problem becomes three-dimensional. Both approximate and 
exact numerical solutions are considered which were applied 
in this case. 

Far-field numerical boundary conditions play an essen- 
tial role in astrophysical applications owing to very large 
length scales usual for these problems. We discuss several 
approaches that may be useful to solve problems similar to 
the stellar wind and interstellar medium interaction. 

1 INTRODUCTION 

The problem of the stellar wind, with the emphasis on the so- 
lar wind, interaction with the interstellar medium has long 
been the topic of interest for astrophysicists and special- 
ists in the field of the solar-terrestrial physics [Q], [g], |pHj, 
[p8[], and [p9|]. In [^|, the application of the continuum equa- 
tions for this problem is systematically discussed. The first 
qualitative model for the interaction of the solar wind (SW) 
and the local interstellar medium (LISM) was proposed by 
E.N. Parker He assumed the interstellar medium to be a 
subsonic stream with the Mach number Moo <§C 1 (see Fig. 1). 

Here HP is a heliopause dividing the S W and the LISM 
flow. Generally speaking the above assumption is not correct, 
since the velocity of the interstellar medium is Voo ~ 20 
km/s and the scattering experiments for the solar radiation 
show that the temperature Too of charged particles consti- 
tuting the LISM is about 10 4 A". Taking into account that the 
number density of charged particles is usually considered to 
be Moo ~ 0.1 cm -3 , the LISM flow can most likely be sup- 
posed supersonic than subsonic. The supersonic model of the 
interaction was first proposed in The important feature 
of this approach lies in the application of the continuum (Eu- 
ler gasdynamic) equations only to the charged particles of 
both counteracting winds. Although the presence of turbulent 
pulsations of plasma is supposed to be insignificant for the 
mean flow structure, their influence is realized by a remark- 
able change of transport coefficients due to the possibility 
of scattering of charged particles on electromagnetic plasma 
fluctuations. This results in the substantial decrease of their 



Figure 1. Schematic picture of the SW-LISM interaction [3] 



Figure 2. General scheme of supersonic interaction 



mean free path compared with that calculated on the basis 
of the Coulomb collisions. The solar wind consists mainly of 
electrons and protons with the number density n e ~ 10 cm -3 
and velocity V, ~ 400-500 km/s and is also considered su- 
personic. The index "e" corresponds to values measured at 
1AU = 1.5 x lO 11 km, that is, at the Earth distance from the 
Sun. Thus, we can consider this problem, from a gasdynamic 
viewpoint, as a an interaction of the supersonic spherically- 
symmetric (or asymmetric) source flow of SW with the uni- 
form supersonic LISM flow. This assumption gave rise to a 
so-called two-shock model. Generally speaking, this model 
can be easily obtained numerically if one choses as initial 
data for this interaction the arbitrary jump between the SW 
and the LISM parameters. 

The source flow is essentially a combination of supersonic 
jet and blunt-body flows which are perfectly well studied 
and described in classical gas dynamics. The extension to the 
solar wind and space situation plj ] , I pBj ] has been validated by 
numerous spacecraft observations of the plasma environment 
of several planets. Different flow regimes and shock- wave 
flow structure for the SW-LISM interaction was discussed 



in[10f 



In \W and [113] the axisymmetric problem of the interac- 
tion was considered on the basis of the shock-fitting approach, 
but due to the limitations of the applied numerical method the 
authors calculated only the upwind part of the flow. In [ |55| ] and 
| p2| ] the problem of the stellar wind interaction with the inter- 
stellar medium was analyzed in the closed region surrou ndin g 
the star. These numerical results confirmed the scheme JlOq], 
but the bullet shape of the internal shock was obtained for a 
broader range of parameters. The general schematic picture 
is shown in Fig. 2. 

Here TS is the inner shock terminating the solar wind within 
the heliopause (HP), BS is the bow shock, or the outer shock. 
TS at a certain point may turn to form the Mach disk (MD). 
At this triple point the reflected shock (RS) and the slip line 
(SL) originate. 

This picture is si milar to that suggested for the SW-comet 
interaction in [106]. The termination shock configuration is 



In ||, [0], [§§], Hi, {jjjl, Jwj ], and [JlO^] the opinion 
was stated that the resonance recharge processes between the 
neutral and the charged particles should strongly influence 
the flow picture. The description of the neutral particle mo- 
tion cannot be made on the basis of the continuum approach. 
For this reason in [^] an approximate method of taking into 
account this influence was suggested and, later, in Jll| ] a self- 
consistent model was developed that takes into account the 
recharge processes. The Monte-Carlo method was used to 
cal culat e the trajectories of neutral particles. As was admitted 
in [104] and confirmed by [p|] and [jllj], the charge-exchange 



caused by its Mach-type reflection from the symmetry axis. 



effect effectively diminishes the Mach number of the LISM 
flow. This justifies development of alternative subsonic mod- 
els of the interaction [p8|. 

Of great importance are also nonstationary problems asso- 
ciated with the variable solar activity. In JTsj], [j75j], and J97| ] 
the problem was investigated of the time-dependent SW per- 
turbation influence on the whole flow structure. The nonsta- 
tionary picture of the flow and the termination shock response 
to the 1 1 -year variation of the solar wind were studied in []17|], 
Q, and [Q. 

T. Matsuda et al. [ p5| | discussed the instabilities of the con- 
tact discontinuity dividing the SW and the LISM flow. These 
instabilities originated in its lateral region and near the stag- 
nation point of the flow. Unstable solutions were obtained in 
the parameter range which was not exactly suitable for the 
solar wind and the local interstellar medium flow and, there- 
fore, were disputed in J96]]. Recently in |^] a hydrodynamic 
instability of the heliopause driven by plasma-neutral charge- 
exchange processes was discussed. 

According to the solar minimum observations by the Uly- 
sses spacecraft, the solar wind properties depend on helioalti- 
tude. The solar wind in this case is no longer spherically- 
symmetric (see also J9^]). In [ |7l| a 3D problem of the in- 
teraction was used and obtained results were compared with 
those for the isotropic solar wind. 

The influence of the interstellar magnetic field can be im- 
portant in the regions, for which the magnetic pressure is 
comparable by its value with the dynamic pressure of the flow. 
Magnetic field causes an increase of the maximum speed of 



small perturbations in the LISM flow, thus leading to the de- 
crease of its effective Mach number Jj3|], [p4|. The ordinary 
Mach number of the LISM is often assumed to be M x — 2. 
At the same time the magnitude of the LISM magnetic field, 
although not known very well, is estimated within 7 x 10 -7 
and 3 x 10 -6 Gauss ft This means, as will be shown later, 
that the magnetic pressure can exceed the value of the ther- 
mal pressure and is the reason of including magnetic field into 
consideration. Its influence can be especially effective when 
the charge-exchange processes are included, since both this 
effects lead to the decrease of the effective Mach number. 
This means that the LISM flow can become sub sonic and the 
bow shock can disappear (see also [^] and [ 104 ]). The stellar 
wind-interstellar medium interaction with taking into account 
magnetic field was first studied in [53], although some of the 
results seem to be misinterpreted (see jl^]). In the latter pa- 
per the problem was investigated by the shock-fitting method 
only in the upwind part of the flow due to the limitations of 
the numerical scheme. In J8(i| ] the solution was presented of 
the axisymmetric problem (the LISM magnetic field strength 
vector was assumed to be parallel to its velocity vector) of the 
SW-LISM interaction in the closed region surrounding the 
star. In [j34j] the shape of the heliopause was studied on the ba- 
sis of the Newtonian approach in the 3D case of the arbitrary 
angle between the LISM magnetic field and velocity vectors. 
The authors found this shape approximately by equating the 
values of the total pressure on the both sides of the heliopause. 
In [ ^T| this problem was first examined numerically. 

The complicated pattern of the flow containing a number of 
interacting shocks, enhanced by various physical phenomena, 
makes it a challenging problem for the application of modern 
numerical methods invented for pure gasdynamic and magne- 
togasdynamic application. Numerical solution of this problem 
is often associated with the solution of such accompanying 
problems as non-reflecting boundary conditions, numerical 
implementation of the condition of the magnetic charge ab- 
sence, etc. These problems are of general importance for the 
young, but quickly developing, field called computational 
fluid dynamics. In Section 2 of this review we present the 
mathematical statement of the problem, write out the system 
of governing equations and boundary conditions. The choice 
of initial conditions for the magnetogasdynamic (MHD) in- 
teraction is discussed. In Section 3 we discuss stationary so- 
lutions of the gasdynamic problem on the basis of shock- 
fitting and shock-capturing methods. In Section 4 different 
approaches are discussed which take into account the charge- 
exchange processes. In Section 5 we describe instabilities 
originating under certain circumstances in this problem and 
the reasons causing them. In Section 6 nonstationary solutions 
are considered resulting from the solar wind disturbances and 
its periodicity. In Section 7 we briefly discuss the effects of 
the solar wind spatial asymmetry. And, finally, Section 8 deals 
with numerical modeling of the solar wind interaction with 
the magnetized interstellar medium. 



2 Mathematical statement of the problem 

To make the paper more concise, we write out in this section 
the mathematical statement of the problem based on the MHD 
equations. The Euler gasdynamic equations can be easily ob- 
tained from the former one by omitting the terms containing 
the magnetic field strength and the equations describing the 
behavior of its components. 

2.1 The system of governing equations 

The system of governing equations for a MHD flow of an 
ideal, infinitely conducting, perfect plasma in the Cartesian 
coordinate system x, y, z, shown in Fig. 2 (v-axis is perpen- 
dicular to the picture plane), can be written as follows (one 
fluid approximation): 
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In system (1) p, u, v, w, B x , B y , and B z are the density and 
the components of the velocity v and of the magnetic field 



strength vector B. We introduced here also the total pressure 
po = p+B 2 /87r (p is the thermal pressure) and the total energy 
per unit volume 
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where 7 = 5/3 is the specific heat ratio corresponding to 
the fully ionized plasma. The quantities of density, pressure, 
velocity, and magnetic field strength are normalized, respec- 
tively, by poo, pooV^o, Voo, and Voo\/Poc, where the index 
"00" marks the values in the uniform LISM flow. Time and 
the linear dimension are respectively related to L/Voo and L, 
where L is equal to 1 AU. The above formulation implies that 
molecular and magnetic viscosities, heat conductivity, and 
anomalous transport effects are neglected. The source term H 
can be both of physical and of geometrical origin and will 
be specified separately for each problem. For example, if the 
flow is axisymmetric system (1) can be rewritten in the plane 
coordinate system x, z as follows: 
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This system is valid in the half-plane xOz and can be ob- 
tained from (1) in the assumption of cylindrical symmetry. Its 
another form is 
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Though both presentations of governing equations are ma- 
thematically equivalent, for numerical reasons it is often more 
convenient to use the former one, since it is supposed to give 
more stable results in the vicinity of the geometrical singular- 
ity x — 0. The Euler gasdynamic equations can be obtained in 
different forms from systems (1) and (2) by assuming B = 0. 



2.2 Initial and boundary conditions 

Calculations are usually performed in the computational re- 
gion between the inner and the outer spherical surface (cir- 
cular in the axisymmetric case). The flow from the Sun is 
supposed to be supersonic at the termination shock distance. 
For this reason we specify all parameter values at the inner 
boundary sphere. The uniform LISM flow is also supersonic 
and we can specify the parameter values at the inflow side of 
the outer boundary. The treatment of outflow boundary can be 
more complicated. In J96| ] all parameters were extrapolated 
with the zeroth order along the fluid particle trajectory. This 
approach cannot be expected suitable for deeply subsonic out- 
flow boundary. Another approach was proposed in |p2||. This 
method is based on introducing imaginary cells next to the 
boundary. These cells are filled with the LISM gas at infin- 
ity. To find the flux through the outer boundary the Riemann 
problem is solved between the imaginary and the adjacent 
cell values. At the boundary segments with the subsonic- 
supersonic transition the rarefaction wave relations are used. 
This results in the following interpretation (see [^] and [|78|]) 
of the method initially developed for a purely gasdynamic 
case and makes possible its extension to MHD problems [|8(i|]. 

Consider the method based on the two strictly nonreflect- 
ing conditions: the well-known extrapolation condition for 
a supersonic exit that provides a characteristically compati- 
ble approximation of the equations on the boundary and the 
procedure J75I, [f78|], developed earlier for gasdynamic flows. 

The idea of application of the relations in the rarefaction 
wave for the realization of the far-field boundary conditions 
lies in the artificial locating of the sonic point on the exit boun- 
dary. If the flow is supersonic at infinity such a procedure gives 
reasonable results and allows one to perform calculations in 
the cases for which other known approaches fail. The interpre- 
tation of our method is the following. Assume that parameters 
inside the chosen computational region fully define the flow 
behavior outside the boundary. In the case of subsonic exit 
the only possible elementary Riemann problem configuration 
for the above system is a rarefaction wave whose fan covers 
the boundary. In this case, if the self-similar variable value 
is known, we can locally continue the internal field to the 
boundary. That is why, an additional condition is that the flow 
velocity attains the sonic value there. 

Consider the hyperbolic system for the vector U of un- 
known variables in the vicinity of the boundary (the right one 



for definiteness) in the form 
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where U = U(x,f), and x is the variable in the direction 
normal to the boundary T, t is time, and A(U) is the coef- 
ficient matrix with a complete set of eigenvectors and only 
real eigenvalues. We seek the solution in the form of a simple 
wave U = U(x,f) = U(£), where £ = - f . By substituting this 
representation into Eq. (4), we obtain 



(A-XI)V ( = 0, A = £, 



(5) 



where / is the identity matrix. Owing to Eq. (5), the vector 
is the eigenvector of A for the eigenvalue A = £. This means 
that we need to solve the following system of ordinary differ- 
ential equations supplemented by the nondifference relation: 



U e = rf(U,A)r(U,A), A(U)=£, 



(6) 



where r is the right eigenvector (the vector-column) of A 
defined up to the scalar multiplier d. The eigenvalue in the 
rarefaction wave varies like A(U) = £. This condition com- 
pletes the system for determining U and d. While realizing 
this boundary condition we must integrate Eq. (6) over £ from 
£o = A(Uo), where Uo represents the initial subsonic param- 
eters inside the region, to £ = £r = 0, that is, to the sonic 
point. Consider this approach, first, for pure gas dynamics. 
Let us choose the vector of unknowns in Eq. (4) in the form 
U = (p,u,v,w,a) T , where p is the density, u is the velocity 
vector component normal to F, v and w are its tangential com- 
ponents, and a is the speed of sound. The minimum eigenvalue 
in this case is A = u — a and the related eigenvector is 
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where 7 is the adiabatic index. System (6) in this case acquires 
the form 
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This system can be exactly integrated, as its invariants are 



0, \u + 



2a 

~ T 



0, v ( = 0, w ( = 0. (9) 



Thus, we obtain 



The index "0" indicates the values belonging to the inner 
region. On the discreet mesh this means that they are taken 
from the center (or from the left side) of the cell adjacent to 
the boundary. Equations (10) must be supplemented by the 
condition at the supersonic exit if (u/a)o > 1: Ur = Uo. 

These conditions are mutually consistent and coincide for 
mo = ao. Note that in this case we did not need the explicit 
expression for d which can be easily found from the second 
and the fifth equations in (8). Besides the exact derivation 
based on the relations in the rarefaction wave, we give here 
the approximate relations keeping in mind such systems for 
which no exact expressions of this kind can be written out. 
For this purpose we first exclude d by substituting the first 
equation from (8) into the other ones. Then we obtain 
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Now, approximating Eqs. (11) by finite differences we ar- 
rive at the following relations: 
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In this approximation solution (12) differs from (10) only 
in the entropy invariant and represents its linearization. 

Now we proceed to the MHD equations. Let us choose the 
unknown vector in Eq. (4) in the form 



U = [p, u, v, w, a, B x , By, B Z ] T . 



(13) 



where, in addition to the purely gasdynamic case, the compo- 
nents appear of the magnetic field strength vector normal (B x ) 
and tangential (B y and B z ) to T. The minimum eigenvalue in 
this case is A = u — a/, where a/ is the largest of the two 
magnetosonic speeds a/ and a s (a/ > a > a s ): 
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The eigenvector corresponding to this eigenvalue is 



7- 1 / 2 
ar = — — r I "0 H -ao I , ur = ar, vr = vo, 



7 + 1 \ 7 - 1 

/ ar\ 
\a J 



w r = wo, pr = po 



(10) 



1, - ^, aB y , aB z , ^ - ^" , 0, (3B y , /3B Z 
P 2p 



(15) 



2pJTTp{a 2 — a 2 s )' p(a 2 — a\ 



In this case system (6) acquires the form 
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Now in system (16) we substitute the equation for a by the 
equation for a/, which can be easily obtained from Eq. (14) by 
direct differencing with respect to £ and by using Eqs. (16): 
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It can be shown in this connection that for any admissible 
values of functions we have i9 + «/ > 0. 

Then, by passing from d to p^, we obtain the reduced 
system of equations that can be approximated similarly to 
Eq.(ll) 
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Note that in contrast to the purely gasdynamic case the 
velocity components tangential to the boundary, generally 
speaking, are different from their internal values in the pres- 
ence of the magnetic field. 

The case of the triple degeneration of eigenvalues (B 2 . + 
Br. — and a 2 — > B 2 /A-np), when $ + af — » 00, can be easily 
avoided by assuming B 2 + B 2 — e, where e is a small positive 
number. 

The approach described above turned out to give stable re- 
sults in contrast to the attempts of a straight applicatio n of the 
well-known non-reflecting boundary conditions, say, [ 102]. It 
is worth mentioning in this connection that a comprehensive 
review (more than 200 references) of various versions and 
modifications of non-reflecting boundary con ditions can be 
found in [^|] and [0| (see also [p7j] and [|io2]]). 



As initial values the jump can be chosen between the SW 
and the LISM parameters at a fixed distance Rf from the 
Sun smaller than the TS stand-off distance. For R < Rf the 
SW parameter distribution is specified. The magnetic field 
pressure in this region is supposed to be negligibly small 
comparing with the SW hydrodynamic pressure, thus B e = 0. 
For R > Rf the uniform distribution of the LISM pressure and 
density is assumed. If we solve an MHD problem, that is, the 
LISM flow is magnetized, it is wise to specify a magnetic field 
strength distribution satisfying the divergence-free condition. 
For this reason the magnetic and the velocity field in the LISM 
and in the SW flow are initially joined in the computational 



region so that B conserves a constant angle with v and div B 
0. This is done by assuming for R > Rf 
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Here (J, V, and W represent the spherical components of 
the velocity vector in the directions R, <f>, and 8, respectively. 
This distribution corresponds to an incompressible fluid flow 
velocity distribution over a sphere, directed along the z-axis 
(#-axis). The magnetic field is initialized in the same way, 
except that the field configuration is rotated about the y-axis 
so that the magnetic field vector is tilted with respect to the 
velocity vector by the desired angle. 

If neutral particles are to be taken into account, we must 
specify their number density and velocity at the inflow. 

3 Stationary gasdynamic calculations 

A gasdynamic model for the SW interaction with the super- 
sonic interstellar wind was first suggested in |Q| . A cylindri- 
cal formulation was adopted which corresponds to a uniform 
LISM flow and a spherically-symmetric SW. The calculations 
were performed on the basis of a simplified thin-layer approx- 
imation. The solar wind was assumed to be decelerated mainly 
in the process of its interaction with the charged component, 
or plasma component, of the LISM. With the latter assump- 
tion calculations can also be made without simplifications 
using the Euler gasdynamic equations. Both shock-fitting and 
shock-capturing approaches can be found in publications and 
we describe them briefly in the following subsections. 

3.1 Shock-fitting methods 

Taking into account the adopted two-shock model with the 
contact discontinuity between the shocks, one can easily use 
shock-fitting methods in the upwind part of the interaction 
region. The idea of this approach lies in the subdivision of 
the computational region into subregions of smooth flow. In 
these smooth subregions any numerical scheme of sufficient 
order of accuracy can be used. Derivatives in this approach 
must never be approximated by finite differences across dis- 
continuities. The latter are traced as boundary lines. Proper 
jump relations are used to determine the change of parameters 
across these boundaries and their new position in the course 
of time. This approach is very economical, since 1) you need 
not perform calculations in the regions of the uniform LISM 
and the spherically-symmetric SW, which are known before- 
hand, thus reducing the size of the computational domain 
and 2) you can avoid spurious oscillations around disconti- 
nuities inherent in application of shock-capturing methods. 
For this reason one can use linear high-order of accuracy 
numerical schemes not worrying about high non-oscillatory 




Figure 3. Self-similar picture of shock wave position [7] 

resolution of discontinuities. In the shock-fitting method their 
position and intensity are determined exactly. Certain diffi- 
culties originate if two or more discontinuities are interacting 
with each other. In this case, in principle, one can still use a 
shock-fitting approach using exact solutions for the problems 
of a discontinuity interaction. The algorithm, however, can 
become rather complicated (see |^], [jsTp, p^l, and |p4]]). 
In the shock-fitting approach a quasi-linear form of the sys- 
tem of governing equations is usually chosen rather than a 
conservation-law form. 

On the basis of the shock-fitting approach the problem 
under consideration was first solved in [7|| by the Babenko- 
Rusanov implicit scheme Due to the limitations of the 
numerical approach only upwind part of the interaction was 
calculated. A polar coordinate system R, 6 was used and cal- 
culations were performed in the region restricted by the ray 
— ^max at which the velocity normal to this boundary 
remained supersonic. The system of linear algebraic equa- 
tions on the computational grid was solved together with 
the Rankine-Hugoniot conservation relations and the rela- 
tions on the contact discontinuity. The computational region 
was located between the termination and the bow shock. The 
steady-state solution was obtained as / — > oo with the boun- 
dary conditions independent of time. 

A physicist usually seeks dimensionless similarity param- 
eters of the problem. For the considered problem they are 
represented by the Mach numbers M e and Moo of the SW 




Figure 4. Discontinuity pattern: shock-fitting approach [113] 



and the LISM flow, respectively, the ratios of their dynamic 



pressures A* = n e V e /n 



and stagnation temperatures 



X — Teo/Tooo- The specific heat ratio 7 is usually adopted to 
be equal to 5/3 which corresponds to the fully ionized plasma. 
If neutral particles are not taken into account, the flow from 
the inner side of TS becomes hypersonic and, therefore, the 
results are only weakly dependent on the choice of M e taken 
at 1 AU. Though results formally depend on \ ( sa Y> f° r X — 1 
the density jump over CD is absent) an analysis of conser- 
vation relations on the discontinuities in the upwind part of 
the interaction region make us conclude that they can be re- 
calculated using results for one particular value of \. As K is 
the similarity parameter (see [^] and for the hypersonic 
solar wind and K > 1 the results are independent of K if the 
distances are measured in units of 1AU x \f~K. In Fig. 3, the 
shapes of the outer (F) and of the inner (G) shock wave, the 
contact surface (S) and some streamlines are presented for the 
case Moo = 2. 



In [ 1 13 1 a new shock-fitting numerical algorithm was de- 
scribed and applied to the problem under consideration that 
provided results with such a high accuracy that in the upwind 
part of the interaction it may represent a testing benchmark 
for all newly developed methods. This method is based on the 
composite explicit-implicit finite-difference scheme [pq]. The 
well-known explicit Lax-Wendroff and the implicit Babenko 
scheme |Q| represent its constituent parts. The algorithm is 
organized in such a way that depending on the CFL num- 
ber value either explicit or implicit approximation is used. 
This approach makes the proposed method more economical 
than a purely implicit scheme. Among the schemes based on 



this approach we can also mention J5^ ] and [j74j]. The latter 
method is the extension of MacCormack's explicit-implicit 
scheme for the steady Euler equations hyperbolic with re- 
spect to one of the space coordinates. A stea dy-st ate solution 
of the SW-LISM interaction was obtained in [ 1 1 3 ] by a quasi- 
marching method in which nonstationary problem was solved 
for t — > oo for each radial ray. Such an approach also allowed 
to save a computational time comparing with a direct solution 
of the nonstationary system. One of the computational results 
is shown in Fig 4. 

One must admit here that, in contrast to the results from 
[Q], calculations were performed for a rather long distance 
in the downwind direction, actually until the flow along the 
z-axis shown in Fig. 4 remained supersonic. This resulted 
in a considerable unjustifiable elongation of the termination 
shock, as will be seen from the subsequent subsection. The 
discrepancy between the results is caused by the limitation of 
the quasi-marching approach that cannot take into account the 
Mach-type reflection of the inner shock from the symmetry 
axis. Nevertheless, the results are without any doubt quite 
reliable up to a certain distance in the wake region. No need 
to mention that the time necessary for obtaining the steady- 
state solution in this case is considerably less than that in the 
case of applying shock-capturing methods. 

3.2 Shock-capturing methods 

In shock-capturing methods we calculate finite differences 
across discontinuities. This may cause spurious oscillations 
of the solution if non-monotone numerical schemes are used. 
On the other hand, all linear schemes of the order of accuracy 
higher than one are non-monotone J38|]. For this reason one 
or another artificial viscosity must be used [|^] or nonlinear 
numerical schemes ought to be applied. It is not our task to 
give a review of high-resolution TVD (tota l var iation dimin- 
ishing) schemes in this paper (see 0] and ||lT|] for a regular 
mathematical background). Although both finite-difference 
and finite- volume methods can be equivalently applied in the 
latter schemes, we shall dwell mainly on the finite- volume for- 
mulation and monotonic upstream schemes for conservation 
laws (MUSCL) approach, since they are more descriptive. 

To solve axisymmetric system (2), let us introduce a polar 
mesh 



Ri=Rt^ + (l-l)AR, 1=1,2,... 
On = (n - 2.5)A0, n = 1,2, . . . ,N; 

AR = (Rmm-Rmm)/(L - 1), 

A8 = tv/(N-4) 



(20) 



with the center in the star. Then for each cell system (2) in the 
finite-volume formulation can be rewritten as follows: 
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(21) 



Here E is the flux normal to the boundary, defined as: 

E = mE + n 2 G, (22) 

where n = (111,112) is a unit outward vector normal to the cell 
surface. 

Equation (21) has a time-discretized conservation-law form 
for an individual computational cell. Various numerical sche- 
mes are specified by the method chosen to calculate the nu- 
merical flux E through the cell boundary surfaces. In J96| ] 
the two-step Lax-Wendroff scheme is used with the second 
order of accuracy. As usual for such schemes, an additional 
smoothing must be introduced to remove high-frequency os- 
cillations and overshoots and undershoots originating near 
smeared shocks owing to the non-monotonicity of the scheme. 
Another possible flux calculation formulas which also include 
artificial viscosity were applied in [f70[]. It is based on the 
ZEUS fractional step code All methods using artificial 
viscosity for oscillation damping contain an empirical vis- 
cosity coefficient which must be adjusted in a way suitable 
for any particular problem. The compromise is between the 
effective monotonization of the solution and its deterioration. 
Nonlinear high-resolution numerical schemes are free from 
this drawback. 

To attain the second order of accuracy in space, a piecewise- 
linear distributio n of parameters inside computational cells 
can be adopted [112]. One can use the simplest "minmod" 
reconstruction procedure 



u: 



1+1/2 



U! 



1 



min mod(AUf +1 / 2 .AUj +3 / 2 ) , (23) 



U/ + i/2 = U ; + - min mod (AIT 



/-l/2>AU; + i/2 



(24) 



minmod(jc,y) = sgn(x) max{0, min[|x|,y sgn(.r)]} , 



where AUj +1 / 2 



u: 



Uf, and Uf +1/2 and Uf +1/2 defined 



by Eqs. (23)-(24) represent parameter values on the right and 
on the left side of the cell surface with the index "/ + 1 /2". 
The index "n" is omitted in these formulas. The reconstruction 
procedure in the angular direction is similar. To attain better 
resolution of the contact discontinuity, one can use more com- 
pressive slope limiting procedure for the density, e. g., 

pf+1/2 = Pi+\ - minmod(Apf +1/ 2,Aft +3/2 ,A) , (25) 



P1+1/2 = Pi + minmod(Api : _ 1/ 2,Apf +1/ 2,A), 



(26) 



(E t „+i/2 + E tB _ 1/2 ) AR + R,ARA6 H,,„ = 0. 



A = 0.25(Ap? +1/ 2 + Ap* + 3 /2 ), 
A = 0.25(Apf_ 1/2 + Apf +1/2 ), 

minmod(x,y,z) = sgn(x) max{0, min[|x|,y sgn(x),|z|]} 

The fluxes E(U R ,U L ) through the cell surfaces can be found 
by different methods. Wide recognition acquired TVD shock- 
capturing methods based on the exact or some of the approxi- 
mate solutions to the Riemann problem. The first steady-state 
solution of the SW-LISM interaction problem in the closed 
region surrounding the star was obtained in |^] on the basis of 
the Osher approximate Riemann problem solver [E4j. In this 



approach an approximate solution to the Riemann problem 
is formed using elementary simple waves which correspond 
to definite eigenvectors of the Euler gasdynamic system and 
separate the regions of the constant flow. In the exact solu- 
tion, a rarefaction wave, a contact discontinuity and a shock 
wave generally appear in various combinations. In the approx- 
imate solution used in [^] shock waves are approximated by 
compression waves. This approach connects parameter val- 
ues on the right and on the left side of the computational cell 
by a number of algebraic relations. Such approach, however, 
seems to be more time-consuming than that based on Roe's 
solution of the linearized Riemann problem | p8| ] (see also 
a comprehensive description of characteristic-based schemes 
for the Euler equations [p^]). A numerical flux in this case is 
calculated as follows: 

E(XJ R ,V L ) = - [E(U L ) + E(U R ) - SIAIS" 1 ^* - U L )] 

(27) 

Here 5(U) and S~ (U) are the matrices formed by the 
right and by the left eigenvectors, respectively, of the frozen 
lacobian matrix 

7 _ dE(U) 

dv 

The value of \J(\J L ,\J R ) is chosen so that the conservation 
relations on shocks are exactly satisfied. The matrix |A| is 
a diagonal matrix consisting of the frozen Jacobian matrix 
eigenvalue moduli. 

The important peculiarity of the latter method is that, al- 
though it gives the solution of the linearized problem, the exact 
satisfaction of the Rankine-Hugoniot relations on shocks pro- 
vides their more adequate and sharp resolution. This method 
was applied to the problem under consideration in p5[ ] and 
@. 

Although all the mentioned methods give essentially simi- 
lar pattern of discontinuities in the computational region, we 
present here, for the convenience of the further discussion, 
the picture from [[7^]. The interstellar plasma number density 
of protons is assumed to be nn+ ~ 1 cm -3 . The velocity of 
LISM relative to the solar system is about 20km/s, while the 
speed of sound of the LISM gas is about 10 km/ s. Thus, LISM 
flow is supersonic. The SW charged particle number density 
is chosen n H+ ~ 10cm~ 3 , its velocity is V e ~ 500km/s, the 
speed of sound c e ~ lOOkm/s at the distance of the Earth's 
orbit (1 AU). The nondimensional parameters of the problem 
are Mach numbers of LISM and SW M x , M e , the relations 
of dynamic pressures K — p e Vf / poo and stagnation tem- 
peratures x = Toe/Tooc of SW and LISM. This corresponds 
to the following values of dimensionless parameters of the 
problem: M x = 2, M, = 5, x = 400, K = 6250. The spe- 
cific heat ratios for SW and LISM are supposed to be 5/3. The 
calculation is performed in the ring region with the inner and 
outer circle radii being R m i n = 10 and R mlix = 500 AU with 
99 and 116 cells in the radial and in the angular direction, re- 
spectively. Constant pressure (below the symmetry axis) and 
constant density natural logarithm contours of the steady-state 
solution are presented in Fig. 5. Results are given in the polar 
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Figure 5. Discontinuity pattern: shock-capturing approach [75] 
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Figure 6. Entropy contours with velocity vectors for the model 
with the Mach number 0.6 [93] 

region with the inner and outer circle radii 10 and 400 AU. 
All peculiarities are seen of the shock wave pattern shown 
schematically in Fig. 2. 

Calculations for a large variety of the stellar wind and in- 
terstellar medium parameters, including those very different 
from the parameters of the solar wind, were performed in J55] ] 
and [ p2[ |. These parametric studies may reflect various situ- 
ations corresponding to ejecting stars and their environment. 
Although from the purely gasdynamic viewpoint the LISM 
flow is supersonic, the presence of the interstellar magnetic 
field and charge-exchange processes can decrease the effec- 
tive Mach number of the uniform interstellar medium flow, 
thus making it subsonic. For this reason several authors [|93|], 
performed calculations of a subsonic interaction on the 
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Figure 7. Entropy contours and velocity vectors for the Mach 
number 5 [93] 
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Figure 8. Geometrical pattern of the interface. Results of the 
numerical calculations for«#oo = 0(1) and Kh oo = 0.14cm — 3 (2); 
curves (3) are the sonic lines. Positions of the bow shock (BS), 
termination shock (TS), heliopause (HP), reflected shock (RS), 
tangential discontinuity (TD), and Mach disk (MD) are shown. 

basis of the Euler equations. It is quite clear that, unlike the 
two-shock Baranov's model, no bow shock can appear in the 
stationary solution. As will be shown later, the termination 
shock in this case does not have a bullet shape and the dif- 
ference between its stand-off distances in the upwind and in 
the downwind direction is not so pronounced. This can also 
be seen in Fig. 6, corresponding to M x — 0.6. If the stel- 
lar wind Mach number is larger than that of the solar wind, 
the termination shock becomes elongated in the backward di- 
rection and the size of the Mach disk diminishes (see Fig 7, 
corresponding to M x = 5). The latter two figures are taken 
from 

Summarizing the subject of this section, we would like to 



admit that, although application of shock-fitting methods can 
provide a solution with a very high precision at low computa- 
tional costs, they can be used only in the flow regions with a 
simple shock-wave pattern which is known beforehand (say, 
the upwind part of the SW-LISM interaction). Attempts to 
promote calculations for larger distances into the downwind 
region re sults in the substantial distortion of the results, like 
it was in Jl 13| ], in which the possibility of the Mach-type re- 
flection of the termination shock was not taken into account. 
A shock-fitting method in its rigorous sense was not real- 
ized, since the exact relations connecting parameters in the 
triple point were not used. The best way out in this case is 
to combine shock-fitting and shock-capturing methods, as it 
was done in jJTT{]. 

4 Stationary solutions including neutral 
particles 

Since we restrict ourselves in this review to the description 
of numerical methods used to solve the SW-LISM interac- 
tion problem, more or less complete description of physical 
processes governing the charge exchange between the plasma 
and the neutral component of the flow lies beyond the scope 
of this work. It can be found in the corresponding references 
mentioned in Introduction. However, it is quite clear that the 
LISM is a partially ionized gas and, therefore, a consistent 
model must be developed accounting for the mutual influence 
of charged and neutral particles. Although both helium and 
hydrogen atoms are present, the authors usually disregard he- 
lium atoms, since their cosmic abundance is much less than 
that of hydrogen atoms. 

The first attempts to estimate the influence of this process, 
say, |0] and |p7||, did not include the influence of the hydro- 
gen atoms on the plasma component. The first self-consistent 
model of the heliospheric interface was suggested in ||s|] . In 
that paper the authors introduced the source terms account- 
ing for the momentum and the energy transfer between the 
two components into the system of the Euler gasdynamic 
equations in the quasi-linear form. These source terms were 
presented in the form: 



pv c {vu - v) 
for the momentum equations (1) and 



p{i - iK 



3kT + 3kT H 
1mn 2mu 



(28) 



(29) 



for the energy equation. Here k is the Boltzmann constant 
and iiih is the mass of the hydrogen atom. For the collision 
frequency between protons and atoms the following formula 
was used: 



Ve = e*zQ, Q = 

m H 



, . 2 , mk(T + t h ) 

v V/y - v) + 



where a is the effective charge-exchange cross-section. 



The hydrogen atoms were assumed to conserve their ve- 
locity and temperature in the process of their interaction with 
protons and the origin of the secondary hydrogen atoms was 
neglected. That allowed to describe the behavior of neutrals 
in the region between the two shocks by the system 

v« = v//^ = const, 
Th = Th^ = const, 
div(p w v//) = -pv c 

Incorporating this source term into implicit method [Q] does 
not represent any difficulty and the solution was obtained for 
various values of dimensionless parameters. The major ef- 
fect of charge exchange on the heliospheric interfaces is to 
decrease the distances to the TS, HP, and BS. The original 
developers of the approach [JsJ] , which considers the flow of 
neutral particles as a hydrodynamic flow, themselves admit- 
ted in ten its main drawbacks: (1) such description is hardly 
justifiable, since the mean free-path of the hydrogen atoms 
is not smaller than the characteristic length of the problem; 
(2) the Maxwellian distribution of the atoms was adopted to 
calculate the plasma momentum and energy losses. This was 
the reason of applying the Monte-Carlo method for simulation 
of the hydrogen atom trajectories It gives a possibility 
to evaluate the source terms in the momentum and the en- 
ergy equations for the plasma component on the basis of the 
kinetic description of the hydrogen atoms. An iteration me- 
thod [0 was proposed to solve both systems of equations. 
The numerical results are discussed in [|ll]] and [fh?]]. This 
approach is nowadays the only one which can be considered 
physically consistent from the viewpoint of computational 
fluid dynamics. In Fig. 8, the geometrical interface between 
the two flows is presented [ pi] ] reflecting the influence of the 
charge exchange processes. The picture is presented in the 
XOZ plane, where OZ coincides with the axis of symmetry 
and is antiparallel to the vector on the LISM velocity (the Sun 
is in the coordinate system origin). The solid and the dotted 
lines in Fig. 8 correspond to n// oc = 0.14 cm -3 and nn x = 0, 
respectively. The parameters of the plasma component were 
the following: n e = 7 cm -3 , V e = 450km/s, M e = 10, 
«oo = 0.07 cm" 3 , Vac, = 25km/s, and Moo = 2. One can 
see a large influence of the described processes on the flow 
pattern and it is quite clear that any realistic calculation of the 
problem must take them i nto c onsideration. 

Recently in [ |7o| | and [115] the simplified approach [^] 
was modified by developing the multifluid description for 
the neutrals which tried to take into account the highly non- 
Maxwellian nature of the neutral distribution. Since in these 
papers the plasma was assumed Maxwellian only in each 
region between the discontinuity surfaces, the neutral popula- 
tion produced by the charge exchange in these regions has the 
same basic characteristics as those of the plasma component. 
Although this approach still remains approximate, one must 
admit that it allowed its developers to obtain the solutions 
of the three-dimensional and nonstationary problems which 
have not been solved yet on the basis of the Monte-Carlo 
method. 
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Figure 9. Entropy contours and velocity vectors [55] (Mao = 15) 

5 Hydrodynamic instabilities in the 
SW-LISM interaction 



As was admitted in [106], the contact discontinuity between 
the inner and the outer subsonic region is potentially the sub- 
ject to the fluid Kelvin-Helmholtz or MHD flute instabilities 
which could produce some of the inhomogeneous structures 
of astropause tails. The study of instabilities accompanying 
the problem of the SW-LISM interaction with and without 
magnetic field is an important independent task which lies 
outside the scope of this review. We are going to clarify here 
only the effect of the choice of numerical methods on their 
origin and development in the problem in which two flows 
possessing substantially different entropies collide at a con- 
tact surface calculated by shock-capturing methods. It is more 
or less clear that discontinuity-fitting methods are hardly ap- 
plicable for this task. 

As was convincingly shown in [ |5^ ] and [63], complicated 
nonstationary patterns can be obtained owing to the Rayleigh- 
Taylor instabilities if the interaction is accompanied by the 
gravitational effect from the star. In the case of the solar wind 
the gravity can be neglected, since the Hoyle-Littleton accre- 
tion radius is about 4.5 AU, which is much smaller than the 
expected location of the termination shock. In [p"5|, however, 
the stellar wind-interstellar medium interaction was studied 
for parameters different from those corresponding to the solar 
wind and the results were obtained clearly indicating insta- 
bility of the flow. The Osher numerical scheme p4j ] was used 
and an algebraically generated irregular grid rather than the 
polar one applied in the previous study of the problem jB3j], 
That type of the grid is supposed to give a finer resolution at 
large distances from the ejecting center. The influence of the 
interstellar medium Mach number was studied in a wide range 
from Mao = 0.6 to 15. A highly nonstationary solution was 
obtained in this case, instabilities originating both near the 
astropause stagnation point and in its lateral region (Fig. 9). 

The Kelvin-Helmholtz instability in the latter region are 
quite possible as the gases with different entropy move along 
the both sides of a contact surface at different speeds. The 
origin of the Kelvin-Helmholtz instability near the stagnation 



point, however, is questionable, since velocities are very small 
there. Such stagnation point instabilities, however, were found 
earlier in the calculations of astrophysical jets J65fa , Js6|l, and 
|p5|]. Similar instabilities were also found in |5f ] and ||57j] 
in the calculations of opposing and forward-facing jets. They 
were also observed in the experiments p6[], [po{, and 
It is worth mentioning, however, that the instabilities of the 
latter type were not obtained in the calculations for parame- 
ters close to the solar wind - interstellar medium interaction, 
as was reported in |p6||. Moreover, N. Pogorelov in his cal- 
culations occasionally obtained instabilities near the stagna- 
tion point, even for a rather rough mesh, using the Roe-type 
MUSCL scheme under unfavorable selection of the parame- 
ter reconstruction method. Note that one or another parameter 
interpolation based on the assumption of the linear, parabolic 
[j^, or higher-order distributions inside computational cells 
is usually used to increase the order of accuracy of the chosen 
numerical scheme. Say, in [ |75| and J76| ] the interpolation of 
characteristic variables was used which had been discovered 
to give much more stable results in the vicinity of the geo- 
metrical singularities at 9 — and 6 — n. Note that both the 
polar grid used in [j76|] and the algebraically generated O-type 
grid incorporated into the algorithm J55| ] possess such a sin- 
gularity. In [ [jc| ] the choice of the form of governing equations 
(2) or (3) was claimed to have a certain effect on the stability 
of results. 

There also exits another aspect of the problem. As was 
mentioned in the previous section, the interaction between 
the interstellar atoms and the solar plasma ions via resonant 
charge-exchange collisions greatly affects the steady-state 
structure of the global heliosphere. The flow of interstellar 
ions is diverted around the nose of the heliopause, whereas 
the neutral particles penetrate into the heliosphere impeded 
aside from the charge-exchange collisions. Thus, there is a ve- 
locity difference between them near the nose which performs 
like a drag force acting on the plasma (see Eq. 28). Since the 
LISM density is larger at the heliopause than that of the SW, 
the contact surface is potentially Rayleigh-Taylor unsta ble. 
That kind of instability was admitted in |p2|] and [115]. Its 



origin was confirmed by several numerical experiments, in- 
cluding the one applying the entirely different particle-in-cell 
code [^fj, and by comparing the instability linear growth rate 
with the theoretical one. However, as was admitted earlier, 
the multifluid hydrodynamic model in which all fliuds were 
supposed to be in a local thermal equilibrium was used in 
the above-mentioned papers to determine the motion of the 
neutral particles. This makes the obtained results questionable 
from the viewpoint of the spatial scale of instabilities. Paper 
| p2| ] also neglects the effect of energetic solar wind neutrals 
created by the charge exchange inside the heliopause. The 
drag force caused by these neutrals tends to compensate the 
inward force described earlier and partially suppress the in- 
stability. 

The study of hydrodynamic instabilities originating in the 
SW-LISM interaction problem is far from its conclusion also 
due to a number of physical phenomena that can affect it, 



such as, the interstellar magnetic field, cosmic rays, etc. It is 
important to admit in this review that, investigating instabili- 
ties, one must be very careful in order to distinguish those of 
physical and of numerical origin. 



6 Nonstationary SW-LISM interaction 

As was mentioned above and is widely accepted, the solar 
wind changes its speed from supersonic to subsonic through 
a termination shock. A number of reasons can cause tempo- 
ral asymmetries of the termination shock and, therefore, the 
interaction pattern as a whole. Among them are disturbances 
of the solar wind and its 11 -year periodicity. Due to these 
reasons the heliospheric shock will move in resp onse to vari- 
ation in upstream solar wind conditions. In [jlOfJ, a kinematic 
analysis is made of the solar wind driven temporal variations 
in the heliospheric termination shock distance. In Jl^] , [jig], 
and Js3| l the motion of this shock was analyzed analytically 
on the basis of one-dimensional gasdynamic model. It was 
admitted that the termination shock would rather be nonsta- 
tionary and would resemble a distorted asymmetric balloon 
with some part moving inward and others moving backward. 
In [Q] the termination shock response to large-scale solar 
wind fluctuations were studied numerically using the Lax- 
Wendroff scheme. The LISM flow was assumed subsonic. 
In ||47|], on the basis of the similar numerical method 11- 
year solar wind variation influence on the inner shock was 
investigated. In J75] ] and [ j76| ] nonstationary problems were 
modeled by a numerical solution of the Euler gasdynamic 
equations (2) in the finite- volume formulation (21) using a 
MUSCL-type TVD high-resolution numerical scheme. The 
suggestion was made of a piecewise-linear distribution of the 
characteristic parameters inside the cells to determine the val- 
ues at their boundaries and slope limiters were used to attain 
a TVD property. The results presented below were obtained 
using the formulas [jj]] : 

E /+1/2 ,„ = E(U L ,U*), (30) 
V L = U,.„ + V',.„AR/2, 
V R = U, + l „ -V' l+u AR/2, 

U;'„ = 5/,„W!,„, 

yyi _ {bm + c)a m + (aj + c)b„, 
a 2 m +b 2 m + 2c 

a = s,;, 1 (u, +1 .„-u,,„), 

b = 5 ; - 1 (U,,„-U ; _ 1 ,„) 

Here c is a small positive value used to avoid division 
by zero. In these formulas S and S _1 are (4 x 4) matrices, 
constructed using right and left eigenvectors of the lacobian 
matrix 9E/9U. 

The fluxes presented by Eq. (30) were defined on the basis 
of Roe's approximate Riemann solver |k3]. 

Fluxes through another pair of cell surfaces can be obtained 
similarly. 



The promotion of the solution in time was performed in 
the following way: 



At dVt 



dt ' 



au (1) 



where t = kAt, k = 0,1, . . . , and At is defined by the time 
resolution and by the CFL condition. 

The LISM proton number density was assumed to be 
iih+ = lcm~ 3 . The velocity of LISM relative to the so- 
lar system is about 20km/s, while the speed of sound of 
the LISM gas is about lOkm/s. The SW protons number 
density was adopted n H + = 10cm -3 , while its velocity is 
V e ~ 500km/s, the speed of sound c e ~ lOOkm/s at the 
distance of the Earth's orbit (1 AU). This corresponds to the 
following values of dimensionless parameters chosen for the 
initial data: = 2, M e = 5, % = 400, K = 6250. The 
stationary initial flow was obtained by the same numerical 
method using a time-stabilization approach. The calculation 
was performed in the ring region with the inner and outer cir- 
cle radii being R m \ n = 10 and i? max = 500 AU. On the inner 
surface all parameters were specified as functions of time by 
formulas: 

(7 = 25, P (t) =p(0)(l+4expHf-3) 2 ]), 
p{t) = 3.2325p 7 (?), 

since this boundary is supersonic. Here U is the radial velocity 
component. Although the choice of parameters is somewhat 
nonrealistic, as far as the solar wind is concerned, it still can 
be used for a qualitative analysis of a nonstationary picture of 
the interaction. 

The initial distributions of pressure (below the symmetry 
axis) and density logarithms are shown in Fig. 5. The size of 
the outer circle in the figures is 400 AU. Similar isolines are 
presented in Figs. 10-13 at different moments of time (in the 
units 1 AU/f/oo). The growing part of the disturbance, inter- 
acting with the inner shock, moves it, first, from the centre, 
Fig. 10, the intensity of IS increasing. Later, it tends to achieve 
the initial position. 

The penetration of the disturbance in the region between 
IS and CD is seen on both charts of isolines at different mo- 
ments of time. The disturbance, while crossing IS, increases 
greatly, and after some time a local pressure maximum orig- 
inates between IS and CD (Fig. 11). This leads to the effect 
of suction of SW gas to the centre accompanying IS motion 
in the same direction. This pressure extremum line gradually 
moves towards CD, see Figs. 12-13. 

The parts of the computational region more remote from 
the source position suffer the same changes later in time. 

As soon as the disturbance meets different points of IS, 
the latter suffers substantial distortion. When the triple point 
on IS is reached, there appear two triple points with dif- 
ferent reflected shocks and recirculation zone between them 
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Figure 10. Pressure and density logarithm isolines [75], t = 8 
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Figure 11. Pressure and density logarithm isolines [75], t = 12 

(Fig. 13). As soon as the source intensity becomes constant, 
the position of IS gradually, but rather slowly, moves to- 
wards its initial position. The maximum increase of MD 
stand-off distance is greater than that of IS along the ray 
9 = 0. 

The relaxation of the flow towards some stationary solution 
is very slow (f > 220). When t > 20, in the region near 9 = 
additional local pressure extrema vanish. The extended vortex 
region originates behind MD and it is up to t = 60, when these 
vortices vanish, probably due to numerical viscosity. 
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Figure 12. Pressure and density logarithm isolines [75], t = 18 
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Figure 13. Pressure and density logarithm isolines [75], t = 21 

Now we present some results fldft concerning the peri- 
odic SW-LISM interaction. In this calculation the LISM pro- 
ton number density, its velocity and the speed of sound are 
«oo = 0.1cm -3 , V x = 20km/s and c x = lOkm/s, re- 
spectively. Parameters of the solar wind change substantially 
within the 1 1 -year period of the solar activity. The concentra- 
tion of charged particles n e = 1.56 cm -3 and the radial SW 
velocity V e = 400km/s are chosen for the minimum of the 
solar activity, while n e — 8 cm -3 and V e = 500km/s corre- 
spond to its maximum (see [™ and [pip). Thus, for dimen- 



rmax= 0. 977 



rm i n= -7. 824 




pmin=-9.210 



Figure 14. Pressure and density logarithm isolines [76], t = 

sionless parameters we have K — 6250, x = 256, — 2, 
and M e — 5 in the minimum and K = 50000, x = 400, 
Moo = 2, and M e = 5 in the maximum of the solar activity. 
These values are chosen as the basic points for the sine func- 
tion approximating the time dependence of K and x within 
1 1 years. The dimensionless time unit is 86.8 days. The time 
step is chosen to be At ~ 0.4 days. Parameter distribution 
corresponding to the stationary solution in the minimum solar 
activity is chosen as initial data. The results are obtained by the 
MUSCL numerical method described above. The calculation 
is performed in the ring region with the inner and the outer 
circle radii being R m i„ — 14 and R mdx — 700, respectively. 

The isolines of the pressure (below the symmetry axis) and 
the density logarithms corresponding to the initial data are 
shown in Fig. 14. They are presented for the outer circle size 
equal to 560 AU. The number of cells is 99 and 116 in the 
radial and in the angular direction, respectively. At the initial 
stage of the flow development, the increase of the parameter 
K results in the TS motion from the Sun. The motion of the 
compression wave through TS causes the origin of a new 
shock wave propagating from the center. This can be seen 
in Fig. 15 corresponding to t = 20. This shock penetrates 
through the contact discontinuity and moves towards the bow 
shock, see Fig. 16 (/ = 36). 

The decreasing part of the periodic function causes a back- 
ward motion of the termination shock towards the Sun. This 
leads to the origin of new flow division surfaces, reverse flow 
zones, and vortices of variable size and intensity in the wake 
region. At t = 60, see Fig. 17, the size of the bullet shape 
termination shock becomes minimum again. Later on a next 
nonstationary shock wave appears, etc. A definite 1 1-year pe- 
riodicity is developed in the shape of the termination shock. 
Parameter distribution between the inner and the bow shock 
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Figure 15. Pressure and density logarithm isolines [76], t = 20 
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Figure 16. Pressure and density logarithm isolines [76], t = 36 

is determined by the propagation of shock waves traveling 
one after another and interacting with the less intensive waves 
reflected from the bow shock (Fig. 18). 

This modeling shows that the flow pattern is substantially 
nonstationary which must be taken into account when analyz- 
ing responses from the space vehicles crossing these disconti- 
nuities. Variable solar activity is usually accompanied by the 
asymmetry of the solar wind which will be considered in the 
next section. 

Application of a high-resolution numerical methods allows 



rmax= 1 . 590 
rm i n= -7. 013 




pmin=-9.210 



Figure 17. Pressure and density logarithm isolines [76], t = 60 
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Figure 18. Pressure and density logarithm isolines [76], t = 84 

one to avoid spurious oscillations near shocks and provides 
their sharper resolution using smaller number of computa- 
tional cells than it is necessary in nonmonotone schemes with 
artificial viscosity. 

7 Nonuniform solar wind - interstellar 
medium interaction 

Distant solar wind deviations from spherical symmetry in- 
duced by the interaction with neutral interstellar hydrogen 



due to photoionization and charge-exchange processes were 
studied in [ pp| ] by a perturbation technique. The penetration of 
neutral particles deep inside the heliosphere results in a sub- 
stantial increase of the distant solar wind temperature. In all 
models described in the previous sections we considered the 
stellar (solar) wind as spherically-symmetric. It can fairly eas- 
ily become asymmetric, since charge-exchange processes are 
clearly more effective in the forward region of the interaction 
because their efficiency is proportional to the velocity differ- 
ence between neutral and charged particles. The perturbation 
analysis has shown that a strongly asymmetric distribution 
of the neutral hydrogen within the heliosphere causes asym- 
metric deceleration and extremely nonuniform distribution of 
the distant solar wind temperature, thus leading to non-radial 
gradients and flows. 

This is, however, only one (external) reason of the solar 
wind asymmetry. Another one is determined by a mere vari- 
ation of the solar wind symmetry in time within its 11 -year 
activity period (see, e.g., speculations in |^]). According to 
the solar minimum observations by Ulysses, the solar wind 
parameters depend on the helioaltitude. These data (see 
and ff^j) indicate that two large polar coronal holes, one in the 
northern and the other in the southern hemisphere, produce 
a hotter, lower-density, higher-speed wind comparing with 
the ecliptic wind. In [ |7l| | the calculations were performed of 
a nonuniform solar wind-interstellar medium interaction us- 
ing the earlier mentioned three-dimensional time-dependent 
ZEUS method [p8|. The solar wind boundary conditions were 
taken from the published Ulysses data. According to them, the 
ram pressure increases 1.5 times from the ecliptic plane to the 
solar pole. That is why, the termination shock in the calcula- 
tions was found to be elongated along the solar axis which, in 
turn, resulted in an increased flow in the ecliptic plane com- 
pared with that over the solar poles. The authors reported a 
pronounced effect of the solar wind asymmetry on the global 
structure of the termination shock and heliopause. Both a two- 
shock (supersonic LISM) and a one-shock (subsonic LISM) 
model were considered. It is worth mentioning once again in 
this connection that, if only the motion of charged particles 
is considered, the LISM flow is definitely supersonic. Pure 
gasdynamic subsonic models for the problem under consid- 
eration are sometimes used to account for the net effect of 
the charge-exchange processes, the influence of cosmic rays 
|^5|], and interstellar magnetic field. The latter subject will be 
discussed in the next section. 

The paper [J7l| definitely indicates, from our viewpoint, 
that any realistic model for the SW-LISM interaction must 
include both neutral particles, the asymmetry of the solar 
wind, and, as a consequence, the solar wind periodicity, since 
the asymmetry varies within the solar activity period. This, of 
course, does not prevent an investigation of different physical 
effects separately. 



8 Solar wind interaction with the magnetized 
interstellar medium 

The presence of the interstellar magnetic field necessitates 
solution of the MHD equations for modeling of the SW- 
LISM interaction. The influence of magnetic field becomes 
important if a magnetic pressure B 2 /87r becomes comparable 
with a dynamic pressure. Magnetic field causes an increase 
of the maximum speed of small perturbations in the LISM 
flow, especially in the direction normal to the direction of 
the magnetic field, thus resul ting in the decrease of the effec- 
tive Mach number J33| ] and rfjjh. Though the magnitude and 
the direction of the interstellar magnetic field is not perfectly 
known, the estimates from [^] indicate the possible impor- 
tance of its presence. The value and the direction of the LISM 
magnetic field can affect the global structure of the interac- 
tion not only directly, but also by modulating the distance 
between the bow shock and the heliopause which determines 
the transparency of this layer for the LISM neutrals. This, in 
turn, affects the observed line profiles of the solar Lyman-a 
backscattered emission and estimates of the LISM parameters 
(see the discussion in [|l|] and J5C|]). If the LISM magnetic 
field is not parallel to its velocity, the problem becomes three- 
dimensional. Later in the section we discuss the results of 
such MHD modeling. 

It is worth mentioning, that although the value of the solar 
magnetic field is often neglected in numerical modeling of 
the problem under consideration at the distances of the termi- 
nation shock, its influence is su bsta ntial at Earth's magneto- 
sphere distances and [ 109]. In [ 10"/ ] the global structure 



of the outer heliosphere was studied in the axisymmetric for- 
mulation for a subsonic interstellar medium with t aking into 
account a time- varying poloidal magnetic field. In [108] and 
Jo7| ] the toroidal magnetic field in the heliosheath was found 
to increase with the distance from the Sun. The complicated 
three-dimensional nostationary behavior of the flow was stud- 
ied which showed the importance of taking this component of 
the magnetic field into account. We restrict ourselves to noting 
this aspect of the problem and will not discuss it below. 

8.1 On the eigenvector and eigenvalue 
systems of MHD equations 

The system of ideal MHD equations in the conservation-law 
form is presented by Eq. (1). This system can be directly 
rewritten in the quasi-linear form which is more convenient for 
the characteristic analysis. It is well-known that the equation 
divB = expresses the absence of magnetic charge. It is 
also evident that, if magnetic charge is absent initially, it 
will not appear mathematically at any time instant. From this 
viewpoint the above equation is excessive. If we rewrite the 
system of MHD equations I ^H 

^ + div p\ = 0, 
at 



97 + (v - 



V)v = -^ 



Vp B x rotB 

p 4np 



dp 
dt 

m 

1st 



dp 
dt 

= rot(v x B) 



in quasilinear form not paying attention to the condition 
div B = 0, the following system is obtained (a is the acoustic 
speed of sound): 
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Solution of the characteristic equation 
det(Ai - XI) = 
gives the following eigenvalues: 

B x 




(31) 



(32) 



(33) 



Note that in one-dimensional treatment the equation for B x 
reduces to 

§ + «f^0, (34) 
at ox 

that is, to the one-dimensional convection equation for B x . 
Of course, in the truly one-dimensional problem (all values 
depend only on the spatial variable x) one can simply as- 
sume B x = const and omit the corresponding equation. On 
the contrary, if we are going to apply the solution of the one- 
dimensional MHD Riemann problem to determine the flux 
through the cell boundary, such an assumption is too exces- 
sive, since only the integral <j> B„da over the whole compu- 
tational cell must be equal to zero. Among eigenvalues (32)- 
(33) the first two correspond to the entropy and B x convection 
waves, A3 4 corresponds to the Alfven, or rotational, waves, 
and the other ones to the slow and to the fast magnetosonic 
wave. Omitting Eq. (34), we reduce the system to 7 x 7. Both 
the extended 8x8 system and the reduced one have real eigen- 
values and a degenerate set of eigenvectors. One can easily 
derive expressions for them. Otherwise, one can refer to |gg], 
and [pit]. As was admitted in fi<| and [0, we can 



use the extended system to derive an approximate solution to 
the MHD Riemann problem. Another possibility is to derive 
formulas for the 7x7 system and use a convection equation to 
find B x (x is normal to the cell boundary) on the cell surface. 
It is also useful to realize that by collecting the source term 
from system (31) to arrive at the conservative form similar 
to (1), we obtain the following system: 



<9U <9E OF dG TI 
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(35) 



where 



H d 



; div 



B (0, 



By By B z v-B 
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This form of the system will be used later to satisfy the 
divergence-free condition. 

Since the application of high-resolution numerical schemes 
to MHD flows and to the problem under consideration, in 
parti cular, has not yet become common (see some tests in 
[103]), we give their description in the next subsection. 



8.2 



High-resolution numerical schemes for 
MHD equation 



TVD upwind and symmetric differencing schemes have re- 
cently become very efficient tool for solving complex multi- 
shocked gasdynamic flows. This is due to their robustness for 
strong shock wave calculations. A general discussion of the 
modern high-resolution shock-capturing methods and their 
application for a varie ty of gasdynamic problems can be 
found in [W2|] and [112]. The extension of these schemes to 



the equations of the ideal magnetohydrodynamics (MHD) is 
not straightforward. First, the exact solution p9| ] of the MHD 
Riemann problem is too multivariant to be used in regular cal- 
culations. Second, several different approximate solvers |E2j], 



@, @], ifjjl, and ||TT4|| applied to MHD equa- 
tions are now at the stage of investigation and comparison. 

The schemes iQ, [||], Q], and [Q, are based on 
the MHD extensions of Roe's linearization procedure In 
|p3|, the attempt of such extension was made and the second 
order upwind scheme was constructed that demonstrated sev- 
eral advantages in comparison with the Lax-Friedrichs, the 
Lax-Wendroff, and the flux-corrected transport scheme |po||. 
Roe's procedure, however, turned out to be realizable only for 
the special case with the specific heat ratio 7 = 2. The reason 
of such behavior of MHD equations is that there is not any 
single averaging procedure to find a frozen Jacobian matrix 
of the system. Another linearization approach is used in p3|], 
fill], [j77j], and [ p4] ] in which the linearized Jacobian matrix 
is not a function of a single averaged set of variables, but 
depends in a complicated way on the variables on the right- 
and on the left-hand side of the computational cell surface. In 
J79| ] and Jj^ this procedure was shown to be nonunique. A 
multiparametric family of linearized MHD approximate Rie- 
mann problem solutions was presented that assured an exact 
satisfaction of the conservation relations on discontinuities. A 
proper choice of parameters is necessary to avoid physically 
inconsistent solutions. 

Consider the one-dimensional system of MHD equations 



dV dF 

dt dx 



0, 



(36) 



where U = yp, pu, pv, pw, e, B y , B Z J and 

/ P u 

pu 2 +po- B 2 x /4n 
puv - B x B y /A-K 
F(U) = puw - B x B z /4n 

e + p Q )u - (uB x + vB y + wB z )B x /4n 
uB y — vB x 
\ uB z — wB x 



In these formulas e = p/(7 — 1) + p{u 2 + v 2 + w 2 )/2 + 
(B 2 + B 2 + S?)/87r is the total energy per unit volume, po = 
p+(B 2 +B 2 +B 2 ) /8-7T is the total pressure,/? and p are pressure 
and density, v = (u,v,w) is the velocity vector, B = (B x ,B y ,B z ) 
is the magnetic field vector, and 7 is the adiabatic index. We 
assume all functions to depend only on time t and on the linear 
coordinate x. Our aim is to construct a solution to Eq. (36) for 
t > for the piecewise-constant initial distribution of U: U = 
Ui for x < and U = Lb for x > 0. It is assumed, owing to 
the divergence-free condition, that B x = B x \ = B X 2 = const. Then 

Let us first find the exact expression for the matrix A — 
AF/AU, where AU = Ui - U 2 and AF = Fi - F 2 . The func- 
tion F being nonlinear, the expression for A is determined 
nonuniquely. In fact, if we choose some nondegenerate sub- 
stitution s = s(U), then from the exact equalities AF = AfAs 
and AU = AuAs it follows that A = Af(Au)~ 1 . The exact an- 
alytic expressions for Af and Ay can be written out explicitly if 
U and F are fractional linear functions of s or polynomials with 
respect to its components. We use here the equivalent trans- 
forms of the type A(BC) = \ {Bx + B 2 )AC + \ (Ci + Ci)AB. where 



The structure and the simplicity of A depends on the choice 
of s. The matrix A is an approximation to the Jacobian matrix 
<9F 

J = — — and must conserve the main hyperbolic properties 
ov 

of J. It must be representable in the form A = Q.rAQ.l, where 
Q.l and Qr are the matrices of its left and right eigenvec- 
tors, respectively, Q.rO,l = /, and / is the identity matrix; 
A = I XjSij\ is the diagonal matrix of the real eigenvalues of 
A, and 5y is the Kronecker delta. 

Then the sought solution to Eq. (36) for t > acquires the 
form 

U(0 = i (U, + U 2 + QgS(C)QtAU), (37) 

F(0 = ^(Fi + F2 + n*|A(£)|£2iAU) (38) 

where £ = f, 5(0 = ||sgn(A, - £)<5,,||, and |A(£)| = 
[[A; sgn(A, — £)<5 t j||. Equations (37) and (38) determine the 
piecewise-constant functions U(£) and F(£) that glue the right 
and the left values of the initial distributions via the system 
of jumps. If the Hugoniot-type condition is valid AF = AAU, 
where A is the jump velocity, then A is one of the eigenvalues 
of A, since AF - AAU = (Af - XAu)As = (A- \E)AuAs 
[A - A£)AU = 0. Thus, det(A - XE) = 0, AU / is the 
eigenvector of A, and relations (37)-(38) describe the jump 
exactly. 

The solution of this kind was first constructed for the equa- 
tions of ideal gas dynamics ||88|]. In [ p2| such solution was 
given for Eq. (36) in the case 7 = 2. The approximate solution 
for an arbitrary adiabatic index was proposed in [[H]]. Here 
we present the procedure for obtaining the extension of Roe's 
linearization procedure for MHD equation and show that it is 
not unique. Thus, the solution from | pd[ ] is the particular case 
of the multiparametric family of approximate solutions to the 
MHD Riemann problem. 

We choose the vector s as a generalization of that for pure 
gas dynamics: 



V 
W 
1 
Y 



( f \ 

y/pW 

\B 



V : 



/R 2 \ 

RU 

RV 

RW 

U s 

RY 
\RZJ 



RT (U 2 + y 2 + W 2 )(j-l) 
Us = 1 

7 27 

_ (Y 2 +Z 2 )R 2 (2-j) 

87T7 

and H — (e + po)/p is the total enthalpy. The vector F in the 
variabless acquires the form 



RU 

F 2 

UV - RYB x /4n 
UW - RZB x /4n 
Ul - {UB X /R +VY+ WZ)B x /4n 
UY - VB X /R 
UZ - WB X /R 



where F2 



(7 - l)i£T/7 - (U 2 + V 2 + W 2 )( 7 - 



l)/2 7 + (Y 2 + Z 2 )R 2 (2 - 7 )/87r 7 . 

Using the new expressions for U and F, we can find the 
matrices Au, A F , and Ax = A F — XAu. We present only the 
expression for A\: 
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where 6i,02,T)i, and 772 are arbitrary parameters. Their origin 
is caused by the presence in the expressions for AF of the 
terms containing the factors ARAY and ARAZ which can be 
attributed both to the terms proportional to AR and Ay or AZ. 
This results in an additional parametrization of the entries of 
the matrices A F and Au. It is not difficult to find that 
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vB x 
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P P 
The following notions are adopted in the above relations: 

P = y/pTpj, U = yfp u/^fp, V = yfpv/y fp, W = y/pw/^P, 

H = ^pH/^p, h y = B y /^p/^p, h z = B z /^p/^p, where 
/ means arithmetic averaging. Besides, 

«=^(^ + ^ + |(Azf 



K = u - A, c = (7 - I) 



u 2 + v 2 + w 2 
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(hy+h 2 z )p 



4irp 4n 
a = 7 (8 - Syhy - S z h z )/2 = 7(5 - q y h y - q z h z ), 

/3 = 7(/i y 5,. + h z 8 z ), 8 = q — — — (hy + h z ) p, 

Sy = 1y ~ l^^A 4 = % - -4^^ 
The equation for a can be rewritten in form 



2-7 
32tt 



ei(AY) 2 + (9i+ V i-l)^^ 
R 



+ - l)('^) 2 + e 2 (AZ) 2 
.ZAZAfi 



,,ZAZAB , ^(ZAR\ 
+ {02 + V2 - 0— 5— + - 1) (-5- J 



The eigenvalues of A are equal to u, u ± £>, where £> = 
iB^I/v^Trp, and to the four roots of the biquadratic equation 



K 4 - 2pK 2 + 2 = 0, (39) 



where 



2p = c 2 + a + b 2 + {h 2 y + h 2 z )p/4n + (3, 
Q = (c 2 + a)b 2 

If c 2 + a > and (h 2 + h 2 z )p/4n + (3 > 0, the roots of 
this equation are real and the diagonal matrix composed of 
the eigenvalues acquires the form 



-, h — v + 2nMsr v , 



A = diag | 



■ a/, u + b, u + a s , u, u — a s , u — b, u — fl/| 



The roots Of and a s are the largest and the least root of 
Eq. (39) ( the fast and the slow magnetosonic waves) and b 
corresponds to the Alfvenic waves. The remaining eigenvalue 
corresponds to the entropy waves. 

The peculiarity of our approach lies in the strict ordering of 
the eigenvalues. This provides the absence of their additional 
nonphysical degeneration which is not inherent in J. Note that 
the choice of other parameter vectors s can break this property. 
In particular, such a degeneration appears if q y and q z are not 
proportional to Y and Z, respectively. This leads to the most 
simple admissible choice of 6: 0\ = 62 = 1. In the MHD 
case, in contrast to pure gas dynamics, it is not possible to 
construct the matrix A depending on a single average vector. 

Let us calculate Q.r and Ql. It is convenient to introduce 
the matrix fl, instead, for which £Ir — AuQ. r . This matrix 
consists of seven columns r. For the eigenvalues \ — u + sa, 
where s = ±1 and a = a/, a s , or 0, the corresponding vector- 
columns r = r(s,a) are the following: 

r = (l, « + 2sa, v — sh y M, w — sh : M, r s , h y N, h z N^j T 
where 

rs = —H + u + v 2 + w 2 + 2sau — (yh y + wh z )sM 
+ [2a 2 -q- (q y h y + qJl z )N] 7/(7 - 1), 



M = M(a) = 



2n(a 2 - b 2 ) 



, N = N(a) = 



a 1 + b- 
a 2 -b 2 



For the eigenvalues A = u+sb(s = ±1) the corresponding 
vector acquires the form 






K 

-K 



\ 



vh* — wh* 
-sh* ^/4n/p sgn B x 
V sh; ^/ 4iv/p sgn B, j 



where ftf = h y /\h\ and h* = /j ; /|h|. 

When using the above formulas for |h| — > 0, the indeter- 
minacies of the type 0/0 must be resolved. This can be done, 
e. g., by the substitution h y = h| sin <^ and h : = |h| cos ip. 

The matrix Q.l can be found similarly by introducing Q./ 
such that CIl = D~ 1 £li, where D is a diagonal matrix specified 
by the equality AuQrD~' Q.i = /. It consists of seven rows. 
For the eigenvalues A = u + sa, where s = ±1 and a — a/, 
a s , or 0, the corresponding vector-rows 1 = l(s,a) are the 
following: 



h = 



7-1 2 
2sn(vr y + wr z )M + (h y r y + h z r z )L, 



h = u- 

7- 1 

k = W + 2nMsr z , Is = — 1, 
k = —r y L + — - , h = -r-L + — 
rv= h L+ 7* 



4n ^(7 — 1) 



,L{a) 



[1+N(a)]p 



4n p(-y - 1)' w 2 
For A = u + sb (s — ±1) we obtain 

Zi = wh* — vh* z , I2 = h = 0, h = h*, U = —hy, 



h = —sh* \J pj4-K sgn 6„ h = s/% \J p/4iv sgn B x 
The matrix D has the form 



D — diag \ \d(af), 2, d(a s ), — d(0), d(a s ), 2, d(ar) 



where 



d{d) — — - 



2 1 2 1 

+ c + a 



[{h], + h\)p/4-K + p](l+N)N/2} 



In practice, if B x — > or |h| — » 0, the indeterminacy 
of the type 0/0 must be resolved in the above relations for 
M — M(a) and N = N(a) at a — a s . This can be done by 
using the biquadratic equation for the roots. It is not difficult 
to find that 



N{a s ) = - 
M{a s ) = - 



(aj -b 2 )(c 2 + a + 2p + b 2 



e\h\ 2 (aj +b 2 ) 



{a} 



& 2 )sgnB. v {c 2 + a)p 



where 



(2-7)(A^) 2 



4-k 



[(2- J?1 )(/<) 2 + (2-r 72 )(/ ; :) 2 ]} 



To resolve the indeterminacy, eigenvectors r(±l, a s ) and 
1(±1, a s ) are multiplied by |h|, and the corresponding d(a s ) 
by jh| 2 . Then, the substitution is made similar to that used for 
the regularization of the Alfvenic eigenvectors. 

No nonphysical degeneration occurs for 7 = 2, since in 
this case a = and /3 = and the characteristic equation has 
only real roots. If we choose 61 = 0% = 1 and 771 = 772 = 2 
for an arbitrary 7, see ||4l|], one can easily show that a > 
and j3 = 0, thus giving only real roots of the characteristic 
equation for any admissible right and left values. 

Another choice is 61 — 82 = 1 and 771 = 772 = 0. In this 
case /3 > 0, whereas a > only for the right and the left val- 
ues connected via the Hugoniot-type conditions. The possible 
degeneration of the formula can be avoided by the regular- 
ization. It is worth mentioning that the above manipulations 
related to the matrix A\ were performed for arbitrary q, q y , 
and q z . This allows us, by substituting q* for q, where 



q = \q — h y q y — h z q z \ + h y q y + h z q z , 

and leaving q y and q z unchanged, to conserve the properties 
of J universally for any < rji — 772 < 2. For small jumps 
this regularization introduces the error of the second order of 
smallness and does not prevent its application for constructing 
the solutions of the second order of accuracy. It does not distort 
the solution for the MHD jumps, since in this case a > 0. 
We also have /3 > for the above choice of jji and 772. Such 
an approach preserves the jump relations and is intermediate 
between the techniques using the exact expression for A — 
AF/AU [||L pjl, and pj] and approaches using different 
approximations to A. 

The family of approximate solutions to the MHD Riemann 
problem presented here generalizes the known approximate 
quasi-linearized and linearized solutions of this problem and 
preserves the Hugoniot-type relations on the jumps. By using 
proper reconstruction techniques, we can increase the order 
of accuracy of obtaining the fluxes in Eq. (38). In this case the 
indices "1" and "2" must be attributed to the parameter values 
on the right and on the left side of the computational cell. 

Due to the complexity of formulas, this approach is much 
more cumbersome than Roe's linearization method in pure gas 
dynamics. Recently in [ p6| ] , a nonlinear approximate Riemann 
problem solver is suggested in which all the waves emanat- 
ing from the initial discontinuity are treated as discontinuous 
jumps. That is why, it is applicable only for weak rarefactions. 
Moreover, the solver proposed is somewhat time-consuming 
and sensitive to the initial approximation for the iteration pro- 
cess. 

Taking into account the above-mentioned remarks some 
simplified approaches are welcome which should (1) satisfy 
TVD property and (2) be enough economical and robust. 

In [[l4j], the second order of accuracy in time and space 
TVD Lax-Friedrichs type scheme is suggested that gives a 
great simplification of the numerical algorithm in the finite- 
volume formulation comparing with the schemes which use 
the precise characteristic splitting of Jacobian matrices. The 
results obtained by this scheme were compared with those 
from J2^ ] and Jz/] ] and a good agreement was observed. In 
this scheme we substitute the diagonal eigenvalue matrix in 
Eq. (38) by the diagonal matrix with the spectral radius (the 
maximum of eigenvalue magnitudes) of A on its diagonal. 
Using the proper parameter reconstruction to find their values 
on the computational cell surfaces, we obtain the second order 
of accuracy. 

This scheme is less dissipative than the original Lax- 
Friedrichs scheme and can be applied to calculations of dis- 
continuous MHD flows. Being much simpler than the scheme 
based on Roe's linearization method, it still gives numerical 
results with the reasonable accuracy. 

Another important subject of discussion is that certain 
initial- and boundary-value problems can be solved non- 
uniquely using different shocks or different combinations of 
shocks, whereas physically one would expect only unique so- 
lutions. The situation differs from that in pure gas dynamics, 



where all entropy-increasing solutions are evolutionary and 
physically admissible. This means that the necessary condi- 
tions of the well-posedness for the linearized problem of their 
interaction with small disturbances are satisfied. In MHD case, 
on the contrary, the condition of the entropy increase is neces- 
sary, but not sufficient. Only slow and fast MHD shocks turned 
out to be evolutionary, while intermediate (or improper slow) 
shocks are to be excluded [ p3| and [p9|], Nonevolutionary 
shocks are not simply unstable in ideal MHD. Their decay 
into evolutionary jumps occurs under infinitesimal perturba- 
tion within infinitesimal time. On the other hand, numerical 
viscosity (including the presence of finite conductivity) and 
numerical dispersion of a numerical scheme make such inter- 
mediate structures existent for a certain time interval before 
their destruction. It is very important to realize that this in- 
terval has nothing to do with the real interval of existence 
of intermediate waves in non-ideal plasma and is grid- and 
numerical scheme-dependent. If the viscosity and/or the con- 
ductivity are substantial, such compound structure can exist 
for a long time depending on the amount of viscosity. This 
was admitted a nd us ed for the explanation of certain physical 
phenomena in [111]. 

In [p^], the ideal MHD Riemann problem was solved with 
initial data consisting of two constant states lying to the right 
and to the left of the centerline of the computational domain. 
Being solved as a strictly coplanar problem, it included a 
nonevolutionary compound shock. Such shocks must decay 
and are not realizable in physical problems. The peculiarity of 
MHD is that there exist discontinuities that are nonevolution- 
ary only with respect to Alfvenic (rotational) disturbances. 
That is why, if a strictly coplanar problem is considered (ve- 
locity and magnetic field vectors lie in the same plane and the 
system of MHD equations includes only two vector compo- 
nents) the construction of the solution is possible both with 
evolutionary and nonevolutionary shock waves. The solution 
in this case is nonunique. If the full set of three-dimensional 
MHD equations is solved and a small tangential disturbance 
is added to the magnetic field vector, a rotational jump splits 
from the compound wave and it degrades into the slow shock. 
This means that the compound wave is unstable against tan- 
gential disturbances and is nonevolutionary in three dimen- 
sions (see [0]). That means that one must be very careful 
reducing the dimension of the system (say, in axisymmetric 
problems) to avoid the origin of nonadmissible solutions. The 
necessity of three-dimensional consideration of MHD Rie- 
mann problems has been admitted recently in One must 
take into account this feature of MHD equations in the con- 
struction of numerical algorithms. In other words, if we are 
going to solve an axisymmetric problem, in general, a full set 
of three-dimensional equations must be used. 

8.3 Shock-fitting calculations 

In [ p^ ] the interaction between the solar wind and the magne- 
tized plasma component of the local interstellar medium were 
calculated by solving the axisymmetric MHD equations using 



the shock-fitting approach. The solar wind was assumed non- 
magnetized. The magnetic field of the interstellar medium 
was assumed parallel to its velocity vector. The system of 
MHD equations in the quasi-linear form was solved by the 
iteration method. Each iteration consisted of the two steps. At 
the first step, the gasdynamic part of the system was solved as- 
suming magnetic field known from the previous time step. At 
the second iteration step, the Maxwell equations were solved 
assuming velocities known. At the first step of each itera- 
tion magnetic field was assumed zero. The g asdy namic part 
was solved by the time-stabilization method [113] (see Sec- 



tion 3). For the known distribution of gasdynamic parameters 
the magnetic field in the considered case can be determined 
analytically as 

pv 



B = B 



J OG ' OC 



These equation follows from the induction equations and 
the MHD jump conservation relations Jjlj ] for B || v. 

The following solar wind and LISM parameters were cho- 
sen: 



= 7 cm 3 , V e = 450km/s, M e = 10, 
= 0.07 cm~ 3 , Voo = 25 km/s, Moo = 2 



The dimensionless value of magnetic field was specified 
via parameter a = 2^/tt/A, where A = Voo/-\/Boo/47rp is 
the Alfven number. 

The geometrical pattern of the flow is shown in Fig. 19. 

It is quite clear from this picture that both the LISM mag- 
netic field and the charge-exchange processes play an impor- 
tant role in the interaction. One can notice that, although a 
shock-fitting approach perfectly well works in the forward 
part of the interaction, its limitations reveal themselves in 
the tail part of the flow. The applicability of the described 
approach is also limited to axisymmetric problems. 

8.4 Shock-capturing calculations 

All mentioned in the previous subsection necessitates the ap- 
plication of shock-capturing methods for modeling of the 
SW-LISM interaction. In |p9j] this problem was solved in 
the closed region surroun ding an ejecting star on the basis of 
the flux-splitting method p6Q. The choice of parameters was 
far from those adopted nowadays for the considered problem. 
Although, as a whole, a physically consistent results were ob- 
tained, the resolution of discontinuities was rather poor and 
some of the obtained data were misinterpreted and disputed 
in&^ 

InjpOJ] the solution of the axisymmetric problem (the LISM 
magnetic field strength vector is assumed to be parallel to its 
velocity vector) is presented for realistic SW and LISM pa- 
rameters. The numerical method applied was developed by 
Pogorelov and is the high-resolution second-order of 
accuracy version of the Lax-Friedrichs scheme. This me- 
thod gives a drastic simplification of the numerical algorithm 
comparing with the methods based on the exact characteris- 
tic splitting of the Jacobian matrices in the MHD equations. 




Figure 19. (1) Geometrical pattern of the interface for a = 2.2 
[13]. Positions of the bow shock (BS), the termination shock (TS), 
heliopause (HP), and sonic line in the solar wind. (2) The same lines 
for a = 0. The discontinuities calculated in [1 1] in the presence of 
the neutral hydrogen are shown by dotted lines 

Simple but very effective numerical boundary conditions at 
the far-field were suggested (see Section 2), which allowed 
one to avoid the influence of spurious reflections from the 
subsonic outer boundary. An effective procedure of satisfying 
the condition of source-free magnetic fields (divergence-free 
condition) was used, which is related to the approach [p3||. 

The formulation of the problem was given is Section 2. 
Following JM]], the numerical flux at the radial cell interface 
(27) was calculated by the formula 

E/+1/2,,, = j [E (Uf + i/2,„) + E (Uf +1 / 2 ,„) + */+l/2.n] , 



(40) 



*;+l/2,n — — R/+1/2.H (uf + i/ 2 ,„ — Uf +1 / 2 ,„) 



Here R; + , /i „ is the diagonal matrix with the same elements 
on its diagonal equal to the spectral radius r (the maximum of 
eigenvalue magnitudes) of the Jacobian matrix |jy : 



4a 2 b 2 R 



r=\U\+a,, 4 = ^({a*f + s/^*y 

b R = B R /(4np)' /2 , B 2 = B 2 R + B 2 e , 
(a*) 2 = (■yp + B 2 /4n)/p, a 2 = 1P /p, 

where if is the radial velocity component and Br and Bg are 



the radial and the angular component of the magnetic field 
strength vector. Similar formulas can be written for the fluxes 
in the angular direction. 

Having the second order of accuracy, the proposed scheme 
is much less dissipative than the original Lax-Friedrichs me- 
thod and provides incomparably better shock resolution. Cal- 
culations were performed in the polar computational region 
with R mm = 24 and R max = 1200. The mesh was R x 9 = 
99 x 116. 

Mathematically, if we choose the initial distribution of pa- 
rameters with divB = 0, the Maxwell equations will pre- 
serve this value in the steady solution. In fact, for problems 
solved numerically the regions of div B / can be accumu- 
lated, especially in the vicinity of discontinuities, see Jl9||. 
It is quite clear that the application of the one-dimensional 
Riemann problem solvers, implying that the component B„ 
of the magnetic field vector normal to the boundary is con- 
stant, contradicts to the condition § B„da — over the whole 
computational cell. In [ fj3| ] an approximate Riemann problem 
solver for MHD equations was proposed on the basis of the 
modified system (35) which is conservative only in a steady 
state. Using the Lax-Friedrichs-type scheme, we do not need 
such modification, since we do not solve the Riemann prob- 
lem to find the fluxes at the cell surfaces. The correction in 
our case can be made by adding to the source term of Eq. (2) 
the value proportional to div B: 



H' = divB 



n fi> B- v B 

0, — , -r^> ~ ; — • "» w 

47T 47T 47T 



This term acts to annihilate the error accumulated if system 
is solved in the conservation-law form (see Jl9[]). It is worth 
noting that this approach lies in the framework of Powell's 
procedure and the correction term is not small only in the 
regions of comparably large errors in div B. This term, in fact, 
is equal by the value and opposite by the sign to the appropriate 
terms proportional to div B appearing by differentiating E and 
G in Eq. (1). The divergence of magnetic field strength vector 
can be approximated over the finite volumes as 

divB = {Rj +l/2 (B„) l+l/z „ + Rl l/2 (B„),_ l/z „]/RjAR + 
[(B„)i, n +i/2 sin6»„ +1/2 + (£;.)/,«- 1/2 sin0„_,/ 2 ]//f/ sin0„A0 

This simple procedure gives a powerful tool for the real- 
ization of the divergence-free condition and necessitates only 
a slight modification of the existing codes for solving MHD 
equations. It seems limited, however, only to steady-state cal- 
culations. 

Numerical results were obtained for the same flow param- 
eters as in the previous subsection and for the three different 
values of the Alfven number: A = 10, 2, y2. The first value 
corresponds to the very small magnetic field (see Fig. 20) and 
can be used for comparison with purely gasdynamic results 
in which the Roe-type Riemann problem solver was used to 
obtain the steady state solution, see Subsection 3.2. In Fig. 20, 
the pressure (below the symmetry axis) and the density loga- 
rithm isolines are presented (A = 10). All features of the flow 
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Figure 20. Pressure (below the symmetry axis) and density 
logarithm isolines, A = 10 
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Figure 21. Pressure (below the symmetry axis) and density 
logarithm isolines, A = 2 

pattern (see Fig. 2) are sharply resolved, while the low time- 
consumption and the simplicity of the algorithm are quite 
clear. In Fig. 21 the same distributions are shown for the case 
A = 2. 

In Figs. 22-23, the case with A = y/l is presented, which 
corresponds to /Jmagn/Pthemai ~ 1.67 at infinity, that is, Boa ~ 
2.5 (dimensional value ~ 2.3 x 10~ 6 Gauss). In Fig. 23, 
the streamlines (lower half) and the magnetic field lines are 
shown. 

Near the symmetry axis the heliocentric distance of the he- 
liopause increases due to the magnetic field tension, whereas 
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Figure 22. Pressure (below the symmetry axis) and density 
logarithm isolines, A = 1.414 
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Figure 23. Streamlines (below the symmetry axis) and magnetic 
field lines, A = 1.414. 

it decreases at the side parts due to the magnetic pressure. 
On the other hand, although the magnetic field is absent in 
the region inside the heliopause, its action is revealed by the 
increased value of the total pressure at infinity. This leads to 
substantial decrease of the termination shock stand-off dis- 
tance in the downstream region. We can see the absence of 
the Mach disk structure in the backward direction. This means 
that the velocity along the termination shock becomes sub- 
sonic. The same occurs if the charge-exchange is taken into 
account JTTj] . The bow shock stand-off distance along the sym- 
metry axis becomes smaller than in the absence of magnetic 



Figure 24. Density logarithm isolines in the symmetry plane 
field. 

As the magnetic field increases, the effective Mach number, 
generally speaking, diminishes. The similar effect is produced 
by the presence of neutral hydrogen atoms and, if both of these 
factors exceed a definite value, the bow shock can disappear. 
In this case system (2) becomes elliptic and other methods 
must be applied for its solution. It is worth noting, however, 
that the speed of propagation of magnetosonic waves depends 
on the direction with respect to the magnetic field vector. 
That means that the effective Mach number varies along the 
bow shock. As will be shown later, this results in a highly 
asymmetric shape of the bow shock if magnetic field is not 
parallel to the velocity vector. 

If the magnetic field strength at infinity is equal to its 
probable upper limit of 3 x 10~ 6 Gauss, the LISM flow be- 
comes subsonic, even if the charge-exchange processes are 
neglected. This can be seen from the simple estimate. We 
can calculate the LISM effective Mach number as Mj? ff = 
PooV^ /(7p 00 + B^/Atx) (this holds, e. g., at the symmetry 
axis for Voo -L Boo • Thus,l/M c 2 ff = 1/M^ + 1 /A 2 and for 
Moo = 4 the effective Mach number remains larger than unity 
only for A > 1. 15. Thus, both effects are of great importance 
for the interpretation of data obtained in the space experi- 
ment. The performance of the algorithm for the realization 
of the magnetic field source-free condition can be seen from 
Fig. 23. The magnetic field lines remain parallel to the stream- 
lines. One can see that no magnetic field penetrates into the 
heliosphere. This result can not be achieved without a special 
treatment of the divergence-free condition. 

The result of the magnetic field influence on the whole 
structure of the flow is not sufficiently investigated. It is clear 
that the flow pattern for Boo If Voo is three-dimensional. 
MHD modeling of the heliopause shape on the basis of the 
Newtonian approximation was performed in |B4J. The results 
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Figure 25. Density logarithm isolines in the plane 
tp = 90° - 270° 

of the three-dimensional modeling of the solar wind interac- 
tion with the magnetized interstellar medium were presented 
in U. 

Calculations were performed in the spherical computa- 
tional region with R mm — 24 and R mM = 1200. The mesh 
isR x 6 x <p = 99 x 116x21. All parameters are the same 
as in the previously described axisymmetric calculation and 
A = 1. In Figs. 24-25 the logarithm density isolines are 
shown for the case of the angle between Voo and Boo equal 
to a = 45° in the cross-sections ip = 0-180° and 90° -270°, 
respectively. The LISM influences in this case the shape of 
the bow shock as well as the location of the stagnation point 
at the heliopause surface in a way consistent with the simpli- 
fied study [34]]. In addition, the results show the existence of 
magnetosheath current layers providing proper rotation of the 
magnetic field with respect to the velocity vector from 45° 
to or 180° at the heliopause surface. The bow shock wave 
stand-off distance is larger in the regions with a larger angle 
between the magnetic field and the shock normal. The contact 
surface is substantially contracted by the magnetic pressure 
in the xy-plane (Fig. 25) rather than in the symmetry xz-plane 
(Fig. 24). As was mentioned earlier, the size of the zone be- 
tween the bow shock and the heliopause is very important in 
view of the charge-exchange processes in this zone. 



interaction. The problem is rather complicated even in gasdy- 
namic and MHD formulations, since the flow pattern contains 
a number of intersecting discontinuities and is substantially 
three-dimensional. The study of this and other astrophysical 
and industrial problems has recently summoned the extension 
of high-resolution numerical methods to magnetohydrody- 
namic flows. The peculiarity of the problem considered in 
this paper is that a continuum approach is applicable only to 
the plasma component of the both flows. Trajectories of neu- 
tral particles must be calculated using a direct Monte-Carlo 
simu lation. Although several authors (see brief discussion in 
[ 1 1C ]) argue that the application of a fluid approximation neu- 



trals gives results close to those obtained on the basis of the 
self-consistent Euler-Boltzmann formulation, there is still a 
necessity to realize this possibility and to perform a critical 
comparison of numerical data for the same set of defining 
parameters. 

The influence of the cosmic rays on the interaction is rather 
well known, but their inclusion into available numerical algo- 
rithms is at the initial stage. 

We examined in this review mainly numerical aspects of 
the problem. This work by no means can be considered as 
an exhaustive description of physical phenomena which take 
place in the interaction region. Discussing the subject, we in- 
evitably had our preferences and paid them more attention. 
The list of references can be extended by hundreds of pub- 
lications on the considered problem, but we hope that even 
those mentioned in this review give an opportunity to realize 
the main processes defining the problem and approaches to 
their numerical modeling. 
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9 Conclusions 

In this work we presented a review of the application of numer- 
ical methods to modeling of the stellar wind interaction with 
the interstellar medium. This is only one among various do- 
mains of their application to space simulation problems. The 
environment of distant stars is not so well investigated as the 
solar system. That is a particular reason of an extremely inten- 
sive study of the solar wind and the local interstellar medium 
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